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ABSTRACT 


The method of general perturbations in rectangular coordinates is 
the most direct of all methods of expansion of the perturbations into 
series, because it is intimately associated with the computation of 
ephemerides. In addition, unlike the method of variation of elliptic ele- 
ments, the method of coordinates does not have the zero eccentricity as 
a singularity. Brouwer T s theory of the general perturbations in the 
rectangular coordinates makes use of the variation of elements in the 
canonical form. However, if the perturbations are being developed into 
trigonometric series with purely numerical coefficients, the use of 
canonical elements is not of any advantage. This fact was recognized 
by Davis who re-wrote Brouwer ! s formulas in terms of the standard 
elliptic elements. Davis T formula contains two terms of order -1 in 
the eccentricity. The presence of these terms causes considerable nu- 
merical inconvenience in the case of nearly circular orbits. We sug- 
gest here use of the Eckert- Brouwer formula for the orbit correction 
as a foundation of a planetary theory. The application of this formula 
leads directly to a vectorial expression for perturbations which is free 
from the disadvantages mentioned above and is also convenient for the 
numerical computations. 

The method of iteration is suggested in computing the effects of 
higher orders. The inclusion of the higher order terms is important 
not only in the planetary case, but also in the case of artificial celestial 
bodies moving in orbits in cislunar space, far away from the Earth. 
Such bodies in their motions resemble more planets or comets than 
satellites. 
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ON SOME POSSIBLE MODIFICATIONS IN BROUWER’S THEORY OF 
THE GENERAL PERTURBATIONS IN RECTANGULAR COORDINATES 


by 

Peter Musen 

Goddard Space Flight Center 


INTRODUCTION 

In our time, the planetary perturbations can be obtained with the same high accuracy by de- 
velopment into trigonometrical series as they can by using step-by-step numerical integration. 

This is the main cause for the revival of interest toward the general perturbations. In addition, 
an expansion of perturbations into series provides us with the possibility of a deeper insight into 
the nature of resonances in the motion of the artificial and of the natural celestial bodies, and 
this is a second reason for further pursuing the analytical or the semi- analytical methods of solu- 
tion in celestial mechanics. The artificial celestial bodies moving in orbits in cislunar space, far 
away from the Earth, resemble a planet or a comet in their motions more than they resemble a 
satellite. In the earth-moon system, the accumulation of the long-period and of the resonance 
effect is greatly accelerated as compared to the speed of accumulation of similar effects in the 
planetary system. Thus the artificial bodies will provide us with an excellent check of our theories 
and will stimulate their further development. Of all methods of expansion of the perturbations into 
series, the method of general perturbations in rectangular coordinates seems to be the most di- 
rect; it is intimately associated with the computation of planetary ephermerides. Unlike the method 
of variation of elliptic osculating elements, this method does not have the zero eccentricity as a 
singularity. The first modern approach to the problem can be found in the work by Brown and Shook 
(Reference 1). A further extension of this idea is due to Brouwer (Reference 2). More recently, 
works using different approaches to the problem were published by Danby (Reference 3) and by 
Musen (References 4 and 5). 

The central idea of Brouwer 1 s theory is one form of variation of astronomical elements. This 
form requires the computation of Langrangian and Poissonian brackets, not for the osculating ele- 
ments, but for the constant elements. Thus this method removes the influence of the variability of 
elements on the coefficients of the disturbing force components from the differential equations, and 
shifts this influence to the modified disturbing force. In the present paper, we express this modi- 
fied disturbing force in terms of Faa De Bruno (Reference 6) differential operators, as we did in 
a planetary theory of a different form (References 4 and 5). The perturbation effects of higher order's 
are transferred from the elements to these operators. The application of the Faa De Bruno oper- 
ators leads to the decomposition of the disturbing force in terms of multipoles with the momenta 
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equal to the perturbations in the position vectors of the planets. In his work Brouwer used the 
system of canonical elements. The formula he deduced resembles the formula for the perturba- 
tions in the coordinates obtained by Brown and Shook (Reference 1). However, if the perturbations 
are being developed into trigonometric series with purely numerical coefficients, then the use of 
canonical elements is not of any advantage. Evidently, this fact was recognized by Davis (Refer- 
ence 7), who transformed Brouwer T s formula and re-wrote it in terms of the standard elliptic 
elements. Davis 1 formula contains, however, two terms of order -1 in the eccentricity. The 
presence of such terms causes considerable numerical inconvenience in the case of nearly cir- 
cular orbits. We suggest here a modification of Brouwer's formula which is free from this disad- 
vantage. In the classical planetary theories the relative positions of the orbital planes are taken 
into account by use of trigonometrical formulas. Such a system is definitely going out of fashion 
because it causes asymmetry in the development of higher order perturbations. Programming for 
digital computers very strongly favors the use of vectors and matrices. For all of these reasons, 
we discard the use of the canonical elements even as an intermediate step, and we suggest, in- 
stead, the use of the Eckert-Brouwer formula for the orbit correction as a foundation of a plane- 
tary theory (Reference 8). This formula in fact represents a solution of the variational equation 
of the two-body problem in terms of position and velocity vectors and also in terms of Gibbsian 
vectors and leads directly to a vectorial expression for perturbations which is both free from the 
disadvantages mentioned above and is convenient for numerical computations. 


BASIC DIFFERENTIAL EQUATIONS 

Here the case of two planets influencing each other is considered. Generalization to the case 
of a planetary system does not present any theoretical difficulty. The differential equation of the 
disturbed motion of the planet with mass m, as disturbed by the planet with the mass m', can be 
written in the form 


d 2 (r + !>r) _ 2 y 

" dt 2 



+ Tc 2 m' V' 




(1) 


where 


M 2 = k 2 (1 + m) . 


By taking into account that 


, M 2 V I (2) 

dt 2 


2 


I 



and by making use of the operators 


D = V 

exp S r * V, 

(3) 

b> 

ii 

Q 

exp & r’ * V' , 3- n d 

(4) 

D" = V' 

exp Sp ■ V # t 

(5) 


we can rewrite Equation 1 as 


= M 2 (D - V) i + k 2 m' f D 'I-D»lV ( 6 ) 

dt 2 r \ r PJ 

By taking the identity 

8r • V V 1 = - JL_ (l-lll 

r r 3 \ r 2 



into account we deduce that 


d 2 hr ^ /j? 3 r r 

"dt*" *1* 


(V 


where 


F = M 2 (D - V - Sr-W) I + k 2 m' ^D' i - D" l) . ^ 

We make use of Brouwer's idea and apply the method of variation of constants to solve Equation 7, 
but avoid the use of canonical variables. The solution of Equation 2 for the undisturbed motion 
has the form 


r = r(t; c r c 2 , c 3 , c 4 , c 5 , , 


where c. are the constant elements. The general solution of the variational equation 


1 

d28r + ^l (l - 3rr '\ • 8 r = 0 
dt 2 r 3 \ r 2 / 

(7a) 

has the form 

6 



, V Sr 

Z_i J Sc 

j-i J 

(9) 


3 



where be. are the arbitrary constants of integration. The method of variation of constants leads 
to the equations 


Sv = 


L 


be. 


j 


3v 

3^T * 

j 


Z 


dSc i 3r 
d t 3c. 

j 


0 


and 


L 

i = i 


dSc i 3v 
dt 3c. 

j 


F • 


It follows from Equations 11 and 12 that 



j=i 



where [c £ , c.] are the Lagrangian brackets 

r -I __ 3r 3y 3r # 3v 

j 3c. 3c. 3c. 3c. 

1 j j 1 

By solving Equation 13 with respect to dSc^dt we obtain 


( 10 ) 


( 11 ) 


( 12 ) 


(13) 


(14) 


d S c. 

i 

d t 


D 

z 

j s l 


(c., c.) ~ ■ F. 
3c. 
j 


(15) 


At first glance. Equations 13 and 15 have the same form as for osculating elements. However, 
the substantial difference between the standard method of variation of the osculating elements and 
Brouwer’s idea is that the Lagrangian brackets in Equation 13 and the Poissonian brackets in 
Equation 15 are formed with the constant elements and not with the osculating elements. The 
sums c. + Sc i are not the osculating elements either. This means that the influence of the per- 
turbations is removed from the elements and is transferred to the operators (Equations 3 through 
5). This transfer produces considerable simplification in the computational procedure, as com- 
pared to the procedure for determining the general perturbations in the osculating elements. 


4 


l■lll■l■lllllll■ n mu mu mi n i 



From Equations 9 and 15 we obtain 


5 r 


£&<>*.>* tt 

i-i 1 1 = 1 j - 1 1 J 1 


Fdt 


(16) 


where (^cj are the constants of integration. In order to shorten the number of trigonometric 
series we use the Hansen and Hill device. We distinguish between the time t under the integrals in 
Equation 16 and between the time contained in the undisturbed position vector r(t) and in its deri- 
vatives which stand outside the integrals. We designate temporarily by r the time associated with 
the undisturbed motion and write r for t(t). The mean anomaly nr +^ 0 associated with r is 
designated by l. We consider r and ? as constants until the integration is performed. After the 
integration, we replace i by t again. Then taking 


(c r c*) - - Oj, c.) 


into account we can write Equation 16 in a symmetrical form 


c 

E 


)r 

3c 


(se. ) 




(t, O* Fdt, 


where r is a dyadic defined by 


I (t, O A 


E E 

, = i , -- 1 


(c , c ) 


dr 

-)c . 


dc. 

} 


3 r 3 r \ 

3 c j 3c J 


In a similar manner we obtain the perturbations in the velocity by, 


(17) 


(18) 


b V 


6 

E 


)v r 


37 


F dt 


(IV) 


On the basis of the principle expressed by Equation 15 each differential equation for the per- 
turbations of an osculating element can be transformed into the equation of our theory simply by 
replacing 30/ 3 C ., the derivative of the disturbing function 0, by 3r/3cj • F. The "disturbing force" 

F is defined by Equation 8. If the element is a vector c, then V c 0 should be replaced by V c r • F. 

In his previous work (Reference 9) the author has obtained the following equations for the osculating 
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P, Q and e: 


dP 

d t 



RR- V Q 0 


/TT? 

na z e 



dQ 

d t 


— 1 — RR • v p n - 

na 2 v'l - e 2 


/TT7 2 

na 2 e 


P 


BO 

Be’ 


de 
d t 


1 - e 2 BO A - e 2 
na 2 e na 2 e 


(P-V Q 0 _ Q-V p 0) . 


From 


r r <f P + r; Q, 


-7 = Cj = a (cos u - e) , and 77 = Sj 


avl-esinu 


we have 


Vp r . I Cj , Vq r = I Sj 


V p r • F - Cj F, and Vq r • F - Sj F, 


and the equations of our theory analogous to Equations 19 through 21 take the form: 


(19) 


( 20 ) 


( 21 ) 


rl<5 P 
cl t 




F - - 


na J v 1 


RR • F, 


( 22 ) 


d.'.g 

d t 


1 - 


dr 

de 


. c t RR • F, and 


1 _ e* 


(23) 


d£ e 

dt 




• F 


(24) 


The coefficients in these equations are formed by using the constant elements. Equations 22 and 
23 can be replaced by the system 


oP dr * P . 

(25) 

&Q . kr > Q, 

(26) 
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and 


d t 


/T 7Z R • F + 1 r R • F 

na 2 e 3e na 2 /l - e 2 


(27) 


A system of approximate equations resembling Equations 25 through 27 appears in Stromgren’s 
theory of special perturbations (Reference 10). In our theory Equations 22 through 24 were exact. 
The theory presented in this paper can also be considered as an extension and improvement of 
Stromgren T s results. 

The classical equations 


da 

_2 an 

and 


dt 

na dl ' 


dl 

1 -e 2 

as) 

2 

dt 

na 2 e 

Je 

na 


do 

da 


have as their analogues in our theory the equations: 


b a 

2 

dl • F, and 

(28) 

a 

na 2 

dl 


1 - e 2 

dr 


(29) 

na 2 e 

le 

2 )a / 

na / 



Each formula for the differential correction of the orbital elements is, in fact, a solution of the 
variational equation (7a), and each such solution can serve as a foundation of a perturbation theory. 
However, the orbit correction formula by Eckert and Brouwer (Reference 8), 


^ r 


tV . r 


dr 


M„ , .21 


a 

a 


+ 


dr 

do 


e , 


(30) 


where 


dr 

a - — 
da 


• in(t 


V 


u 


(31) 


— = H r + K- 3 ! 


(32) 


cos u + e 
1 - e 2- 


(33) 


7 



and 


2 


(34) 


K = 


1 - 



sin u , 


is associated in a most intimate way with the determination of the general perturbations in r and v. 
Direct application of this formula, by-passing the use of the canonical elements, leads immediately 
to a form of r which is in agreement with the general spirit of this theory. 

Substituting Equations 24 and 27 through 29 into Equation 30 we obtain: 



The expressions 

rxRitl-e 2 — and r * R + v'l - e 2 

contain the eccentricity as a factor which will be canceled with the eccentricity in the denomina- 
tor. Our form of P makes it evident that e = 0 is not a singularity. The use of the standard el- 
liptic elements and of Gibbsian vectors leads to a numerical theory which is valid also for nearly 
circular orbits. If the numerical theory of the general perturbations were to be developed for the 
osculating elliptic elements directly, then the eccentricity would appear in the denominator in 
several places, thus causing considerable difficulty in the case of nearly circular orbits. 

Substituting 


c 1 P + Sj Q, Cj a(cos u - e), s x ~ a v 1 ~ 


e" sinu , 


— = c P + s 0 Q, c„ - - - sin u, and s ~ + — /\ - e 2 
H 2 2 2 r 2 r 


(36) 


into Equation 35 we obtain the following decomposition of T: 


r p p 

u 


r„ P0 + r 91 OP + r 99 0 0 


r 33 RR 


(37) 
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where 


vi .. e 2 


<*4 C 3 - ^3 


2 


( C e C o - C 0 C ) , 


Vl - e 2 _ 2 _ 

" + - 2 — (C 4 S 3 ~ S 4 C 3> Cc 5 S 2 -C 2 S ), 

na 2 na 2 


r 21 = + 


v 1 - e 2 


(s 4 c 3 s 3 c 4 ) + (s q c, - s 0 cj 


5 v 2 3 2 5 ' 


r _ v 1 - e 2 

22 ' + J 3 " a 3 a 4 „ 

n a* na* 


C ® 4 s 3 - S 3 S 4 ) + (s 5 s 2 - s 2 s 5 ), 


33 ^ C 1 S 1 - S 1 °A’ 

na 2 vl - e 2 


c 3 • - A. il_ i 1 _ e 2 sin 2u, s 3 a * — (cos u - e) cos u 


c 4 - - a (1 + f. s i n 2 u ) , 


l' 1 - c 4 


• — ( cos u - e) sin u, 


3 . . a 

cos u - e + — n (t - t Q ) _ sin u 


. and 


a v 1 - o 2 


3 >. a 

sin u - — . n (t - t Q ) — cos u 


The developments of c. , s. into the trigonometric series in * with numerical coefficients can 
be obtained by using either the analytical methods or numerical Fourier analysis. The series for 
c., I. are obtained from the corresponding series for c i? s t by replacing i by f. After we ob- 
tained the developments of c i , s i? c\, s., then the developments of F into a double Fourier series 
with the arguments l and i can be performed on an electronic computer without any great difficulty. 


Setting 


P F _ F v Q * F = F 2 , and R * F - F 3 , 
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( 38 ) 


dWj 

~dt 


r f 

hi 


+ ^12 F 2 


dW 0 


dT - r 21 F 1 + r 22 F 2 - and 


(39) 


dWa 

d t 


r f 

33 3 


(40) 


we deduce 

= (W 2 P + W 2 Q + W 3 R) r=t + (br) 

= P bg + 0 *?? + R , (41) 

where 


(^r) - dT x r + — b J L n + a — 
dl 0 Ba 


ba 

a 


Br 

Be 


be , 


and the symbols 


bY 


• 


ba 

— > and be 
a 


now designate the constants of integration. The expressions for the partial derivatives are given 
by Equations 31, 32, and 36. 


DECOMPOSITION OF THE DISTURBING FORCE 

The system Equations 38 through 41 can be solved either by developing br into a power series 
in the masses m and m' or by the method of iteration. If we choose the first way we set: 

Sr = bj r + S 2 r + S 3 r + * • • • , 

br' = r' ■+ b 2 r' + b 3 r' + • • * • , 
bp - b 2 p + b 2 p + b 3 p 4 . . . . t an( j 

b k p = \r' - \r, 
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where S k r, S k r' and S k p are the kth order in m and m' . We have 


exp (Sr * V) = T 0 + T a + ? 2 + T 3 + t 

exp(Sr' -V') =T' + TJ + T' + T' + . . . . , and 
exp(Sp -V') = T" + T" + T" + T" + ...._ 

where T k , T k ' , and T k " (k = 0, 1 , 2 , . . .) are Faa De Bruno differential operators. By setting 


8 k = S k r.V, 


we have 


T 0 = 1 , Tj = 8,. T 2 = s 2 + i 8 2 . and T 3 = 8 3 + 8, S 2 + - Sf 


By replacing S k in Equation 42 by 


S' = S u r' * V\ 

k k 


or by 


k = \sp.v', 


(42) 


we obtain the operators T k ' and T k ", respectively. 

In order to obtain the expansions of the operators 

T, , T\ T" D, D' , and D" 

k k k ’ 

in powers of perturbations one can use formulas for the potentials of multipoles. By setting 

p = r - a 

and designating the multipoles moments by a n a 2 , a 3 , . . . , and the del-operator with respect to 
r by v, we have: 


a 


i 



a i *P 
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n n 1 3 1 

a i * v a 2 * V “= + — p • a i P ' a 2 “ — a i * a 2 ’ 

p p 5 p 3 


ij* Va 2 • Va 3 - = 


1 

P 


15 3 / 

" — P * a ! P ■ a 2 P ‘ a 3 + — (p ■ a j a 2 * a 3 + p ‘ a 2 a 3 ' a i + P ‘ a 3 a i ' a 2> ' 

P P 5 

a 1 * V a 2 • V a 3 * V a 4 • V 1 - 


( 43 ) 


105 15 , 

+ — P ’ a ! P • a 2 p ’ a 3 p * a 4 " — (p • a 2 p • a 3 3j • a 4 + p • a 3 p ■ a, a 2 • a 4 

P P 7 


+ P ' 3 ! P ’ a 2 a 3 ‘ a 4 + P ’ a j P ‘ a 4 a 2 • a 3 + P • a 2 p ‘ a 4 3j 


a 3 4 p • a 3 p • a 4 a i • a 2 ) 


^ <«1 * a 4 

p 5 


They remain valid if one of the moments is replaced by a matrix. Equations 38 through 41 re- 
quire a knowledge of the components of F in the directions of P,Q and R. In order to obtain these 
components we have to replace one of the moments by P, 0 or R respectively. 

In a planetary theory we have to replace 

a k \ r ' \ r '’ an d * k P 


and 


V can mean either V or V' . 


The application of the relations 43 to the expansion of Equation 8 was fully discussed in the 
author's previous work and therefore the details can be omitted in the present article. On occa- 
sion it might be more convenient to use the process of iteration instead of expanding br and br 
in powers of masses. Replacing one of the moments in Equation 43 by the idemfactor and setting 
the remaining moments equal to *>r, we have 


6r ■ V V I : i A r r' h - -J_ b r , 


12 



1 ( S r • V) 2 vi = — (r - Sr) 2 r + — (r • Sr Sr + i r Sr 2 ), 


IS 

r 2? 


i (Sr- V) 3 Vi = + ( r - Sr) 3 

6 r 2r 9 


15 

77 


[r r* Sr Sr 2 + (r* Sr) 2 Sr] ; 


From the last equations and from Equation 3 we deduce that 


where 


(D - V - 8r ■ V V) A r + BS r, 
r 


A = - il a 2 + JL/3 2 

2 r 7 2 r 5 


35 

+ CL - - 


15 


2 r 7 


a'/ 3' 2 + 


R 3 15 2 

Jo = + — a a + 


= r * Sr , and P 2 - Sr 2 . 


In an analogous way we deduce that 


D' -L = A' r' + B' Sr', 


1 3 15 ~'2 . 3 ^'2 , 35 a' 3 - 15 a/3' 2 + 


A' = — r— + a’ a.' - + 

2 r ' 5 2 r 


2 r 


2 r 


B' 


1 3 # 15 t2 

+ a a * + 

r' 3 r' 5 2 r* 7 


a' = r' . Sr\ 


ft 2 = Sr' 2 , 


D" i = A" p + B" Sp , 


(44) 


(45) 

(46) 


(47) 

(48) 


(49) 


(50) 
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-H a" 2 +_L S'' 2 + 

35 «3 15 „ 2 

a s a /3 1 + 

2p~ 

2p s 

2 P 9 2p 7 

1 + 3 a" 15 

« 2 

CL «$**••• 

p 3 p s 

2p 7 



a" = p • 8p , 211^ 


( 51 ) 


/3 " 2 = bp 2 . ( 52 ) 

The expansions of A and B into double series in a and fi 2 converge very fast. The negative powers 
of r, r' and p appearing in the coefficients can be expanded into double Fourier series, with the 
arguments t and V , by means of harmonic analysis. By making use of Equations 44, 47 and 50 
we can write Equation 8 in a compact form, 

F=(ArfBi)r)+k 2 m' [(A' r' + B' b r' ) + (A" p + B" Sp) ] , ( 53 ) 

which is convenient for computing the perturbations by iteration. Decomposing Equation 53 along 
P, Q and R we obtain the expressions for F lf F 2 , F 3 to be used in association with Equations 38 
through 41. In the author r s earlier theories (References 4 and 5) Sr was decomposed along r, v 
and R or along r, Rxr, R. The same form of F (Equation 53) can also be used conveniently in these 
earlier theories if the method of iteration is preferred. 


CONCLUSION 

The results given in this article represent a modification and extension of the results by 
Brouwer, Davis and Stromgren. Equation 8, the basic integral equation, which determines the per- 
turbation Sr, can be solved by means of iteration with the help of Equations 35 and 44 through 53. 
The solution by expanding the perturbations into power series with respect to the masses can also 
be achieved using the system of operators given in the author's previous work. We suggest the 
application of the double harmonic analysis to expand the negative powers of the mutual distance. 
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Appendix 

Notations 


k - the Gaussian constant 

m - the mass of the disturbed planet; the mass of the Sun is one 

fj? - k 2 (1 + m) 

m' - the mass of the disturbing planet 
r - the undisturbed position vector of the planet m 


u 

V 

P, Q, R 


n = /j,a 


* 3/2 


4 - n t + 


- the undisturbed eccentric anomaly of m 

- the undisturbed velocity of m 

- the Gibbsian vectors of m 

- the undisturbed semi- major axis of m 

- the undisturbed eccentricity of m 

- the undisturbed anomalistic mean motion of m 

- the mean anomaly at the epoch of m 

- the undisturbed mean anomaly of m 

- the undisturbed position vector of m' 


r + hr ~ disturbed position vector of m 

S k r - the perturbations of the kth order in the position vector of m 
Sr - the total perturbations in the position vector of m 
Sv = dS p/dt - the perturbations in the velocity of m 

r' + Sr' - the disturbed (actual) position vector of m ' 

8 k r' - the perturbations of the kth order in the position vector of m' 
Sr 1 - the total perturbations in the position vector of m' 
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r - r 


P = 

P = 

Sp = Sr' - Sr 

I - the iden factor 

V - the del-operator with respect to r 
V' - the del-operator with respect to r' 

D" = V' exp Sp *V' 

D' = V 7 exp Sr' * V 7 
D = V exp Sr * V 

V p - the del- operator with respect to p 
V Q - the del- operator with respect to 0 
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